## Replication Code for:

## Immigration Policies and Access to the Justice System:
## The Effect of Enforcement Escalations on Undocumented
## Immigrants and Their Communities 
## Authors: Reva Dhingra, Mitchell Kilborn, Olivia Woldemikael






####Load required packages and functions####

rm(list=ls())
# install required packages if they are missing:
list.of.packages <- c("tidyverse", "stargazer","scales",
                      "lfe","xtable", "MatchIt", "choroplethr",
                      "lubridate")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])];
if(length(new.packages)){install.packages(new.packages)}
lapply(list.of.packages, require, character.only = TRUE)

## Package version check
packinfo <- as.data.frame(installed.packages(fields = c("Package", "Version")))
packinfo[which(packinfo$Package %in% list.of.packages),c("Package", "Version")]

# choroplethr choroplethr   3.7.0
# lfe                 lfe 2.8-5.1
# lubridate     lubridate   1.7.9
# MatchIt         MatchIt   3.0.2
# scales           scales   1.1.1
# stargazer     stargazer   5.2.2
# tidyverse     tidyverse   1.3.0
# xtable           xtable   1.8-4

sessionInfo()
# R version 4.0.2 (2020-06-22)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Catalina 10.15.6
# October 17, 2020


## Load functions
'%!in%' <- function(x,y)!('%in%'(x,y))

## Set working directory to whatever directory CRIME.R is in
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

####Load and prepare necessary data####

## Load county-level crime dataset merged with 
## county level demographic variables.

## Crime data taken from Uniform Crime Reporting Data Series
## compiled by National Archive of Criminal Justice Data
## https://www.icpsr.umich.edu/web/NACJD/series/57?start=0&sort=TITLE_SORT%20asc&SERIESQ=57&ARCHIVE=NACJD&PUBLISH_STATUS=PUBLISHED&rows=50
## file:"county_ucr_offenses_known_1960_2017_rda/county_ucr_offenses_known_yearly_1960_2017.rda"
 
## Variables from UCR include ViolentCrimeIndex, PropertyCrimeIndex, TotalCrimeIndex,
## coverage_indicator. All crime indices are in per 100,000

## Demographic data taken from 2012-2017 ACS 5 Year Estimates
## Data downloaded from socialexplorer.com.
## Variables from ACS include BlackPopPct, WhitePct, AsianPct, Male15_34PopPct,
## SingleMotherHHs, UnemploymentRate, ManufacturingIndustryPct, 
## RetailHospitalityIndustryPct, RenterPopulationPct, MedianHomeValue,
## PctHousingUnitsNoVehicle

ucr_dat<-readRDS("ucr_dat2.RDS")


## Add county-level foreign born data

## Data downloaded from National Historical Geographic Information Systems (NHGIS)
## at Integrated Public Use Microdata Series (IPUMS).
## Year:             2013-2017
## Geographic level: County (by State)
## Dataset:          2017 American Community Survey: 5-Year Data [2013-2017, Tracts & Larger Areas]
##    NHGIS code:    2013_2017_ACS5b
## NHGIS ID:      ds234
## https://data2.nhgis.org/main
fb<-read.csv("nhgis0011_ds234_20175_2017_county.csv", stringsAsFactors = FALSE)%>%
  mutate(ForeignBornPct=100*AH8XE013/AH8XE001,## Convert foreign born prop to %
         ForeignBornNonCitizenPct=100* AH8XE021/AH8XE001, ## Convert foreign born non citizen prop to %
         ## Convert foreign born non-citizen from LA to %
         ForeignBornNonCitizenLA=100*(AH8XE026+AH8XE023+AH8XE024+AH8XE025)/AH8XE001,
         ## Extract FIPS code from GISJOIN variable
         fips_state_county=paste0(substr(GISJOIN, 2,3),substr(GISJOIN,5,7)))%>%
  select(fips_state_county, ForeignBornPct,ForeignBornNonCitizenPct,ForeignBornNonCitizenLA)
## Merge foreign born variables to ucr_dat by fips code
ucr_dat<-left_join(ucr_dat, fb, by=c("fips_state_county"))



## Read in ILRC 2017 data
## Downloaded March 28, 2020 using internet archive to access 2017 version
# of the map. Current version of the map is available:
# https://www.ilrc.org/local-enforcement-map
ilrcFrame<-readRDS("ILRC2017Score.RDS")%>%
  mutate(ILRC2017Score=as.numeric(ILRC2017Score),## Convert score to numeric
         ## Create indicator for whether county received score of 5 or higher from ILRC
         ## as sanctuary county.
         SanctuaryCounty=ifelse(ILRC2017Score>=5,1,0),
         SanctuaryCounty=ifelse(is.na(ILRC2017Score), NA, SanctuaryCounty),
         ## Recode ILRC score so that 0-7 corresponds to high to low 
         ## cooperation with ICE
         ILRC2017Score=as.numeric(recode(ILRC2017Score, 
                                         "0"="7","1"="6","2"="5",
                                         "3"="4","4"="3","5"="2",
                                         "6"="1","7"="0")),
         ## Fix dropped leading zeroes.
         FIPS=ifelse(nchar(FIPS)<5, paste0("0", FIPS), FIPS))%>%
  distinct(FIPS,.keep_all = TRUE)

## Merge ILRC data to ucr_dat
ucr_dat<-left_join(ucr_dat, ilrcFrame, by=c("fips_state_county"="FIPS")) 



## Read in and prepare National Crime Victimization Survey
## data from 2013-2017.
## Data downloaded in April 2019 from 
## Bureau of Justice Statistics
## https://www.bjs.gov/developer/ncvs/index.cfm
## Note that new years of data are added annually, so 
## the PERSONAL_VICTIMIZATION dataset will have a
## different number of years available depending
## on when one downloads the data.
ncvs<-read.csv("NCVS_PERSONAL_VICTIMIZATION_2013-2017.csv")


ncvs_dat<-ncvs%>%
  ## Create indicator for whether crime victim identifies as Hispanic
  mutate(Hispanic=ifelse(hispanic==1, 1, 0),
  ## Create indicator for whether crime victim notified police of victimization
         Notify=ifelse(notify==1,1,0),
  ## Create indicator for whether response is in 2016 or 2017
         Year2017=ifelse(year>=2017, 1, 0),
  ## Create indicator for whether crime victim identifies as Black.
         Black=ifelse(race1r==2,1,0))%>%
  ## subset to 2016-2017 data.
  filter(year>=2016)


## Load Secure Communities Data
## To assemble this dataset, we combine county-level
## economic and demographic indicators from the
## American Community Survey with county-level monthly 
## crime data from the UCR and deportation intensity
## data from the Transactional Records Access Clearinghouse
## at Syracuse University.

## ACS variables include Hispanic Pop Pct, Black Pop Pct,
## Median Household Income Logged, Foreign Born Pop Pct

## UCR: Jacob Kaplan has concatenated
## UCR data into a county by month dataset available here:
## https://www.openicpsr.org/openicpsr/project/100707/version/V11/view
## In some counties, crimes are only reported quarterly, biennially, 
## or annually. We imputed monthly values for these counties by averaging 
## the total annual crime across each respective time period, 
## following the approach taken by Miles and Cox (2014). 
## We also acknowledge the caveat made by Maltz and Targonski (2002) 
## who report missing data issues with the UCR stemming from reporting 
## variation within counties.
## UCR variables include ViolentCrimeIndex, PropertyCrimeIndex, TotalCrimeIndex,

## TRAC: In August 2019 we scraped the TRAC website for 
## data on monthly county-level deportations under the Secure Communities
## program. Link: https://trac.syr.edu/phptools/immigration/secure/
## We then normalized the deportation rate to reflect deportations
## per 100,000 people.

SC_Analysis<-readRDS("SC_Analysis.RDS")


#### Plot Figure 1####

## Total Crime Index by County in 2017
data("df_pop_county")## Load county fips map data for list of fips

## Create tibble with 2017 Total Crime Index from ucr_dat
mapplot<-ucr_dat%>%filter(year==2017)%>%
  select(fips_state_county, TotalCrimeIndex)%>%
  mutate(fips_state_county=as.numeric(fips_state_county))

## merge df_pop_county to mapplot for complete list of fips
df_pop_county<-left_join(df_pop_county, mapplot, by=c("region"="fips_state_county"))
df_pop_county<-df_pop_county %>%
  select(region, TotalCrimeIndex)%>%
  mutate(TotalCrimeIndex=ifelse(TotalCrimeIndex>8000, 8000, TotalCrimeIndex), ## Top code crime index at 8000 for ease of presentation
         TotalCrimeIndex=ifelse(is.na(TotalCrimeIndex), 0, TotalCrimeIndex))
colnames(df_pop_county)<-c("region", "value")## Rename TotalCrimeIndex to value for choropleth package


data("continental_us_states")## Load data to set zoom to continental US

## Plot using choroplethr 
county_choropleth(df_pop_county,
                  title =  "",
                  legend = "Total Crimes Per 100,000", num_colors = 1,
                  state_zoom = continental_us_states)

#### Plot Figure 2####
## 2017 Hispanic share of county population
data("df_pop_county")## Load county fips map data

## Create tibble with FIPS and Hispanic Pop Pct
map2plot<-ucr_dat%>%filter(year==2017)%>%
  select(fips_state_county, HispanicPct)%>%
  mutate(fips_state_county=as.numeric(fips_state_county))

## merge df_pop_county to mapplot for complete list of fips
df_pop_county<-left_join(df_pop_county, map2plot, by=c("region"="fips_state_county"))
df_pop_county<-df_pop_county %>%
  select(region, HispanicPct)%>%
  mutate(HispanicPct=ifelse(is.na(HispanicPct)==TRUE, 0, HispanicPct))## Set counties without ACS entry to 0


colnames(df_pop_county)<-c("region", "value")## Rename TotalCrimeIndex to value for choropleth package
data("continental_us_states")## Load data to set zoom to continental US

## Plot using choroplethr 
county_choropleth(df_pop_county,
                  title =  "",
                  legend = "Hispanic Population Pct", num_colors = 1,
                  state_zoom = continental_us_states)


####Plot Figure 3####
## Get 2017 ILRC Scores from ucr_dat
ilrcPlot<-ucr_dat%>%filter(year==2017)%>%
  mutate(region=as.numeric(fips_state_county), 
         value=as.character(ILRC2017Score))%>%
  select(region, value)%>%filter(!is.na(value))

## Plot
choro = CountyChoropleth$new(ilrcPlot)
choro$set_zoom(continental_us_states)
choro$legend = guide_colorbar(ticks.linewidth = .5, title="Enforcement\nIntensity")
choro$render()


####Table 1 Results####

## Total Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total<-felm(TotalCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year>=2016)
# summary(lsdv_1617_total)

## Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop<-felm(PropertyCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016 )
# summary(lsdv_1617_prop)

## Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol<-felm(ViolentCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016)
# summary(lsdv_1617_viol)


## Produce latex output using stargazer package
stargazer(lsdv_1617_total, lsdv_1617_prop, lsdv_1617_viol, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Hispanic Pct X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:LSDV1", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level")

## Substantive effect size interpretation
lsdv_1617_total$coefficients[1]/ucr_dat%>%
  filter(year>=2016)%>%
  summarize(SD=sd(TotalCrimeIndex, na.rm=TRUE))%>%pull()


#### Table 2 Results####

## Total Crime Index: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total<-felm(TotalCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year>=2016 )
#summary(lsdv_1617_total)

## Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop<-felm(PropertyCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016)
#summary(lsdv_1617_prop)

## Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol<-felm(ViolentCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016)
#summary(lsdv_1617_viol)


## Format model output for latex using stargazer package
stargazer(lsdv_1617_total, lsdv_1617_prop, lsdv_1617_viol, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Enforcement Intensity X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:ILRC", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level")


#### Table 3 Results####

## Regress indicator for whether respondent notified police on 
## interaction of Hispanic indicator and 2017 indicator. Use
## NCVS population weights.
## Note that we use this kludgey method to
## extract robust standard errors from felm object,
## which otherwise cannot be passed directly to stargazer.
mod_basic<-felm(Notify~Hispanic*Year2017,
             data=ncvs_dat, weights = ncvs_dat$weight)
mod_basic_stats<-summary(mod_basic, robust=TRUE)
mod_basic_RSE<-mod_basic_stats$coefficients[,2]
mod_basic_PVal<-mod_basic_stats$coefficients[,4]


## Add age, gender, race
mod_demcovs<-felm(Notify~Hispanic*Year2017+factor(ager)+factor(gender)+factor(Black), 
                        data=ncvs_dat, weights = ncvs_dat$weight)
mod_demcovs_stats<-summary(mod_demcovs, robust=TRUE)
mod_demcovs_RSE<-mod_demcovs_stats$coefficients[,2]
mod_demcovs_PVal<-mod_demcovs_stats$coefficients[,4]

## Add city size, urban vs. rural vs. suburban region of US
mod_demcovs_area<-felm(Notify~Hispanic*Year2017+factor(ager)+factor(gender)+factor(Black)+
                                 factor(popsize)+factor(msa)+factor(region),
                          data=ncvs_dat, weights = ncvs_dat$weight)
mod_demcovs_area_stats<-summary(mod_demcovs_area, robust=TRUE)
mod_demcovs_area_RSE<-mod_demcovs_area_stats$coefficients[,2]
mod_demcovs_area_PVal<-mod_demcovs_area_stats$coefficients[,4]


## Produce latex output using stargazer package
stargazer(mod_basic, mod_demcovs, mod_demcovs_area,
          se=list(mod_basic_RSE, mod_demcovs_RSE, mod_demcovs_area_RSE),
          p =list(mod_basic_PVal, mod_demcovs_PVal, mod_demcovs_area_PVal),
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Hispanic Pct", "Year 2017", "Hispanic X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          omit=c("ager",  "Black", "popsize", "region", "msa", "gender", "Constant"),
          column.labels   = c("Probability of Notifying Police"),
          column.separate = c(3),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes","Yes"),
                           c("Demographic Controls", "No", "Yes","Yes"),
                           c("Geographic Controls", "No", "No","Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:NCVS", 
          title = "Difference-in-Difference Linear Probability Model: Change in Probability of Reporting Crime to Police 2016-2017")



####Plot Figure 4####
SC_Analysis%>%group_by(Date)%>%
  summarize(PctCountiesAnyDeportation=sum(Activated,na.rm = TRUE)/n())%>%
  filter(Date>"2008-01-01")%>%
ggplot(aes(x=Date, y=PctCountiesAnyDeportation))+
  geom_line(color="red", size=2)+
  theme_classic()+
  xlab("Year")+
  scale_y_continuous(labels=percent)+
  ylab("Percentage of Counties With At \nLeast 1 Secure Communities Deportation")+
  theme(axis.text=element_text(size=24),
        axis.title=element_text(size=24,face="bold"),
        axis.text.x = element_text(angle = 90, hjust = 1))+
  scale_x_date(date_breaks = "1 year", date_labels = "%Y")


#### Table 4 Results####

## DV Total Crime Index ~ Deportation Rate + County Covs+ County, Month, Year FE, SEs Clustered by County
TOTAL<-felm(TotalCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct|fips+factor(Month)+factor(year)|0|fips,
            data=SC_Analysis)
#summary(TOTAL)

## DV Violent Crime Index ~ Deportation Rate + County Covs+ County, Month, Year FE, SEs Clustered by County
VIOLENT<-felm(ViolentCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct|fips+as.factor(Month)+factor(year)|0|fips,
              data=SC_Analysis)
#summary(VIOLENT)

## DV Property Crime Index ~ Deportation Rate + County Covs + County, Month, Year FE, SEs Clustered by County
PROPERTY<-felm(PropertyCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct|fips+as.factor(Month)+factor(year)|0|fips,
               data=SC_Analysis)
#summary(PROPERTY)

## DV Total Crime Index ~ Deportation Rate + County Covs+ County, Month, Year FE, SEs Clustered by County + Linear Time Trend
TOTAL_T<-felm(TotalCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct+Time|fips+factor(Month)|0|fips,
              data=SC_Analysis)
#summary(TOTAL_T)

## DV Violent Crime Index ~ Deportation Rate + County Covs+ County, Month, Year FE, SEs Clustered by County + Linear Time Trend
VIOLENT_T<-felm(ViolentCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct+Time|fips+as.factor(Month)|0|fips,
                data=SC_Analysis)
#summary(VIOLENT_T)

## DV Property Crime Index ~ Deportation Rate + County Covs + County, Month, Year FE, SEs Clustered by County + Linear Time Trend
PROPERTY_T<-felm(PropertyCrimeIndex~DeportationRate+HispanicPopPct+BlackPopPct+logMedianHHIncome+ForeignBornPct+Time|fips+as.factor(Month)|0|fips,
                 data=SC_Analysis)
#summary(PROPERTY_T)

## Format model output using stargazer package
stargazer(TOTAL,PROPERTY, VIOLENT,TOTAL_T,PROPERTY_T, VIOLENT_T,
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Deportation Rate"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          omit=c("HispanicPopPct", "BlackPopPct", "logMedianHHIncome", "ForeignBornPct", "LogPop", "Time"),
          table.placement = "H",
          column.labels   = c("Total", "Property", "Violent","Total", "Property", "Violent"),
          column.separate = c(1,1,1,1,1,1),
          add.lines = list(c("Month Fixed Effects", "\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                           c("County Fixed Effects", "\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark"),
                           c("Year Fixed Effects", "\\checkmark","\\checkmark","\\checkmark","","",""),
                           c("Linear Time Trend","","","","\\checkmark","\\checkmark","\\checkmark")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:LSDVICE", 
          title = "DiD Estimates: Change in Crime Indices at US County Level 2002-2014 Due to Rollout of Secure Communities")

####SI-A1 ILRC Definitions####

## This section uses xtable to print the table included 
## in SI-A1 detailing the ILRC definitions of county cooperation
## with ICE
tableDescription<-read.csv("ILRCDefinitions.csv")%>%select(1,2)
colnames(tableDescription)<-c("Policy", "Coverage")
print(xtable(tableDescription,
             align=c('cp{2in}p{2in}'),
             label="tab:ILRCDef",
             caption="ILRC Coding: Local Policies on Assistance in Immigration Enforcement",
             floating=TRUE, table.placement="H", type="latex",
             booktabs=TRUE), include.rownames=FALSE)



####SI-A2 Pre-Treatment Trends Table 1####

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2015 and 2016
lsdv_1516_total<-felm(TotalCrimeIndex~HispanicPct:Year2016|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year %in%c(2015, 2016))
#summary(lsdv_1516_total)

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2014 and 2015
lsdv_1415_total<-felm(TotalCrimeIndex~HispanicPct:Year2015|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year %in%c(2014, 2015))
#summary(lsdv_1415_total)

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2013 and 2014
lsdv_1314_total<-felm(TotalCrimeIndex~HispanicPct:Year2014|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year %in%c(2013, 2014))
#summary(lsdv_1314_total)

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2012 and 2013
lsdv_1213_total<-felm(TotalCrimeIndex~HispanicPct:Year2013|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year %in%c(2012, 2013))
#summary(lsdv_1213_total)

## Format model output using stargazer package
stargazer(lsdv_1213_total,lsdv_1314_total,lsdv_1415_total,lsdv_1516_total,
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Hispanic Pct X Year 2013", "Hispanic Pct X Year 2014", "Hispanic Pct X Year 2015", "Hispanic Pct X Year 2016"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("2012-2013", "2013-2014", "2014-2015","2015-2016"),
          column.separate = c(1,1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:LSDV2", 
          title = "Difference-in-Difference Estimates: Change in Total Crime Indices at US County Level 2012-2016")


####SI-A3 Data Quality Sensitivity Analysis Table 1####

## Repeat Table 1 Analyses while subsetting to counties which have a coverage
## indicator of 100. 

## Total Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total100<-felm(TotalCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                         subset = year>=2016 & coverage_indicator==100)
#summary(lsdv_1617_total)

## Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop100<-felm(PropertyCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                        subset = year>=2016 & coverage_indicator==100)
#summary(lsdv_1617_prop100)

##  Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol100<-felm(ViolentCrimeIndex~HispanicPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                        subset = year>=2016 & coverage_indicator==100)
#summary(lsdv_1617_viol100)


## Format output for latex using stargazer package
stargazer(lsdv_1617_total100, lsdv_1617_prop100, lsdv_1617_viol100, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Hispanic Pct X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:Coverage", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level: 100\\% Coverage Counties")


####SI-A4 Foreign Born IV Table 1####
## Total Crime Index: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total<-felm(TotalCrimeIndex~ForeignBornNonCitizenPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year>=2016)
#summary(lsdv_1617_total)

## Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop<-felm(PropertyCrimeIndex~ForeignBornNonCitizenPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016)
#summary(lsdv_1617_prop)
## Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol<-felm(ViolentCrimeIndex~ForeignBornNonCitizenPct:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016)
#summary(lsdv_1617_viol)


## Format model output for latex using stargazer package
stargazer(lsdv_1617_total, lsdv_1617_prop, lsdv_1617_viol, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Foreign Born Non-Citizen Pct X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:FB", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level")


####SI-A5 Top Quintile Hispanic Pop Pct Table 2####

## Get top quintile value for Hispanic Pop Pct in 2016
TopQuintile<-ucr_dat%>%
            filter(year>=2016)%>%
            pull(HispanicPct)%>%quantile(.8, na.rm=TRUE)

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total<-felm(TotalCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year>=2016 & HispanicPct>=TopQuintile)
#summary(lsdv_1617_total)

##Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop<-felm(PropertyCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016&  HispanicPct>=TopQuintile)
#summary(lsdv_1617_prop)

##Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol<-felm(ViolentCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016& HispanicPct>=TopQuintile)
#summary(lsdv_1617_viol)


## Format model output using stargazer package
stargazer(lsdv_1617_total, lsdv_1617_prop, lsdv_1617_viol, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Enforcement Intensity X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:ILRC", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level")


####SI-A6 Matching Analysis Table 2####

## Create 2017 dataset with variables to be matched on
## Only include observations with complete data availability
matchData<-ucr_dat%>%filter(year==2017)%>%
  select(SanctuaryCounty,state_abb,TotalCrimeIndex, ViolentCrimeIndex, PropertyCrimeIndex,
         HispanicPct, CountyPopDens,Male15_34PopPct, SingleMotherHHs,BlackPopPct,
         UnemploymentRate,ManufacturingIndustryPct,RetailHospitalityIndustryPct,RenterPopulationPct,
         MedianHomeValue,PctHousingUnitsNoVehicle,CountyPopulation_ACS)%>%drop_na()

## Perform nearest neighbor matching using MatchIt package
m.out1 <- matchit(SanctuaryCounty ~ HispanicPct+CountyPopDens+Male15_34PopPct+SingleMotherHHs+BlackPopPct+
                    UnemploymentRate+ManufacturingIndustryPct+RetailHospitalityIndustryPct+RenterPopulationPct+
                    log(MedianHomeValue)+PctHousingUnitsNoVehicle+log(CountyPopulation_ACS),
                  data=matchData, method = "nearest", distance = "logit",
                  replace=TRUE)

## Get balance table from matchit output
balanceTable<-summary(m.out1)
## Select relevant columns. Having trouble placing sideways table
## in latex document, reduce number of columns so table will
## fit vertically rather than horizontally
balanceTable$sum.all<-balanceTable$sum.all%>%select(c(1,2,4))
balanceTable$sum.matched<-balanceTable$sum.matched%>%select(c(1,2,4))

## Format names for nice printing
row.names(balanceTable$sum.all)<-c("Distance", "Hispanic Pct", "County Pop. Density",
                                   "Male 15-34 Pop. Pct", "Single Mother Household Pct",
                                   "Black Pop. Pct", "Unemployment Rate", "Manufacturing Industry Pct",
                                   "Retail Hospitality Industry Pct", "Renter Population Pct",
                                   "Log Median Home Value", "Pct Housing Units No Vehicle",
                                   "Log Population")
row.names(balanceTable$sum.matched)<-row.names(balanceTable$sum.all)

## Print balance table for unmatched data
balanceXTableAll<-xtable(balanceTable$sum.all, label = "tab:BalAll", 
                         caption = "Unmatched Data: Treatment and Control Balance")
print(balanceXTableAll)

## Print balance table for matched data
balanceXTablePruned<-xtable(balanceTable$sum.matched, label = "tab:BalPruned", 
                            caption = "Matched Data: Treatment and Control Balance")
print(balanceXTablePruned)


## Select paired data and regress
matches<-matchData[c(row.names(m.out1$match.matrix),m.out1$match.matrix ),]

## Total Crime Index DV: Sanctuary City and controls, SE clustered on state
totalCrimeMatched<-felm(TotalCrimeIndex~SanctuaryCounty+HispanicPct+CountyPopDens+Male15_34PopPct+SingleMotherHHs+BlackPopPct+
                          UnemploymentRate+ManufacturingIndustryPct+RetailHospitalityIndustryPct+RenterPopulationPct+
                          log(MedianHomeValue)+PctHousingUnitsNoVehicle+log(CountyPopulation_ACS)|0|0|state_abb, 
                        data=matches)
#summary(totalCrimeMatched)
## Property Crime Index DV: Sanctuary City and controls, SE clustered on state
propertyCrimeMatched<-felm(PropertyCrimeIndex~SanctuaryCounty+HispanicPct+CountyPopDens+Male15_34PopPct+SingleMotherHHs+BlackPopPct+
                             UnemploymentRate+ManufacturingIndustryPct+RetailHospitalityIndustryPct+RenterPopulationPct+
                             log(MedianHomeValue)+PctHousingUnitsNoVehicle+log(CountyPopulation_ACS)|0|0|state_abb, 
                           data=matches)
#summary(propertyCrimeMatched)

## Violent Crime Index DV: Sanctuary City and controls, SE clustered on state
violentCrimeMatched<-felm(ViolentCrimeIndex~SanctuaryCounty+HispanicPct+CountyPopDens+Male15_34PopPct+SingleMotherHHs+BlackPopPct+
                            UnemploymentRate+ManufacturingIndustryPct+RetailHospitalityIndustryPct+RenterPopulationPct+
                            log(MedianHomeValue)+PctHousingUnitsNoVehicle+log(CountyPopulation_ACS)|0|0|state_abb, 
                          data=matches)
summary(violentCrimeMatched)


## Format model output using stargazer package
stargazer(totalCrimeMatched, propertyCrimeMatched, violentCrimeMatched,
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Sanctuary County"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          omit=c("CountyPopDens","Male15_34PopPct", "SingleMotherHHs","BlackPopPct",
                 "UnemploymentRate","ManufacturingIndustryPct","RetailHospitalityIndustryPct",
                 "RenterPopulationPct","MedianHomeValue","PctHousingUnitsNoVehicle",
                 "CountyPopulation_ACS","Constant", "HispanicPct"),
          column.labels   = c("2017 Total Crime", "2017 Property Crime","2017 Violent Crime"),
          column.separate = c(1,1,1),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on state in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:MatchedSample", 
          title = "Nearest Neighbor Matched Sample: Sanctuary County Crime Rates in 2017")


####SI-A7 Data Quality Sensitivity Analysis Table 2####


## Repeat Table 3 analyses while subsetting to counties which have a coverage
## indicator of 100. 

##Total Crime Index: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_total<-felm(TotalCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                      subset = year>=2016 & coverage_indicator==100  )
#summary(lsdv_1617_total)

##Property Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_prop<-felm(PropertyCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016& coverage_indicator==100)
#summary(lsdv_1617_prop)

##Violent Crime: County and year fixed effects, SEs clustered on County for years 2016 and 2017
lsdv_1617_viol<-felm(ViolentCrimeIndex~ILRC2017Score:Year2017|fips_state_county+year|0|fips_state_county, data=ucr_dat,
                     subset = year>=2016& coverage_indicator==100)
#summary(lsdv_1617_viol)


## Format model output using stargazer package
stargazer(lsdv_1617_total, lsdv_1617_prop, lsdv_1617_viol, 
          no.space = TRUE, digits=3, digits.extra = 0,
          covariate.labels=c("Intensity X Year 2017"), 
          star.cutoffs = c(0.05, 0.01, 0.001), 
          multicolumn = FALSE, 
          model.numbers = FALSE,
          table.placement = "H",
          column.labels   = c("Total Crime Index", "Property Crime Index", "Violent Crime Index"),
          column.separate = c(1,1,1),
          add.lines = list(c("Year Fixed Effects", "Yes", "Yes", "Yes"),
                           c("County Fixed Effects", "Yes", "Yes", "Yes")),
          model.names = FALSE,style="ajps",
          notes="SEs clustered on county in parentheses",
          dep.var.caption = "", omit.stat = c("BIC", "f", "ser", "adj.rsq"), 
          dep.var.labels.include =FALSE,
          label = "tab:ILRCCoverage", 
          title = "Difference-in-Difference Estimates: Change in Crime Indices Between 2016 and 2017 at US County Level: 100\\% Coverage Counties")

